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Abstract 

Suppose  we  are  given  a  video  of  a  rotating  object  and  suppose  we  want  to 
determine  the  rate  of  rotation  solely  from  the  video  itself  and  its  known  frame  rate. 
In  this  thesis,  we  present  a  new  mathematical  operator  called  the  Geometric  Sum 
Transform  (GST)  that  can  help  one  determine  the  angular  frequency  of  the  object  in 
question.  The  GST  is  a  generalization  of  the  discrete  Fourier  transform  (DFT)  and 
as  such,  the  two  transforms  have  much  in  common.  However,  whereas  the  DFT  is 
applied  to  a  sequence  of  scalars,  the  GST  can  be  applied  to  a  sequence  of  vectors. 
Most  importantly,  we  show  that  the  GST,  like  the  DFT,  can  (1)  be  used  to  estimate 
frequency  and  (2)  can  be  computed  surprisingly  quickly.  Indeed,  we  provide  a  Fast 
Geometric  Sum  Transform  (FGST)  algorithm  that  computes  the  GST  in  0(N  log  N) 
matrix-vector  multiplications,  where  N  is  the  number  of  images  in  the  video  sequence. 
This  is  a  vast  improvement  over  the  0(N2)  such  multiplications  required  for  a  direct 
computation  of  the  GST.  The  remainder  of  this  thesis  is  devoted  to  proving  other 
properties  of  the  GST  and  giving  proof-of-concept  numerical  examples. 
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Determining  angular  frequency  from  video 

WITH  A  GENERALIZED  FAST  FOURIER  TRANSFORM 

I.  Introduction 

Consider  a  video  sequence  of  image  frames  in  which  an  object  moves  throughout 
the  sequence,  such  as  those  presented  in  Figure  1.  This  figure  is  an  image  sequence 
of  a  rotating  jet  of  plasma  which  is  being  ejected  through  a  small  propulsive  device 
called  a  thruster.  Suppose  you  want  to  know  the  rate  in  which  this  jet  is  rotating,  but 
only  the  video  sequence  and  its  sampling  rate  is  known.  In  this  thesis  we  introduce 
an  operator  called  the  Geometric  Sum  Transform  (GST)  that,  when  applied  to  such  a 
video  sequence,  allows  us  to  estimate  this  rate  of  rotation.  Specifically,  this  operator 
is  useful  in  finding  the  rate  of  translation  or  alternatively,  the  rate  of  rotation  of  an 
object  captured  on  video.  It  turns  out  that  the  GST  is  closely  related  to  the  Discrete 
Fourier  Transform  (DFT).  As  such,  we  begin  here  by  recalling  the  basic  properties  of 
the  DFT. 

The  DFT  is  an  operator  which  is  widely  used  in  signal  and  image  processing. 
Letting  N  E  N,  the  DFT  is  defined  over  the  following  sequence  space  of  all  discrete¬ 
time,  complex-valued  N-periodic  signals: 

£( Ztv,  C)  :=  {x  :  Z  — y  C  |  x[n  +  N]  =  x[n\,Vn  6  Z}.  (1) 

Specifically,  for  any  x  e  £(ZN,C),  define  the  DFT  of  x  to  be  DFT  (a;)  G  C), 

N—l 

DFT  (x)[n\:=J2e~2-^x[n'}.  (2) 

n'= 0 

There  are  three  main  reasons  why  the  DFT  is  popular:  (1)  it  determines  fre¬ 
quency,  (2)  there  is  a  fast  algorithm  to  compute  it,  and  (3)  it  helps  us  understand 
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Figure  1:  This  15- frame  image  sequence  shows  a  rotating  jet  of  plasma  coming  out 
of  a  200- Watt  Hall  thruster,  which  is  a  small  device  intended  to  propel  a 
satellite.  Not  only  is  the  plasma  ejected  from  the  thruster  rotating,  the 
brightness  of  the  plasma  is  also  not  constant.  These  images  were  collected 
by  Liu  as  discussed  in  [8],  but  the  true  rate  of  rotation  of  the  thruster  was 
an  unknown  parameter.  Applying  the  GST  to  a  sufficiently  well-sampled 
video  sequence  will  provide  a  good  estimate  of  the  true  rate  of  the  plasma’s 
rotation. 

filters.  The  GST  generalizes  the  first  two  of  these  properties  to  vector-valued  func¬ 
tions.  We  introduce  this  new  operator  by  generalizing  the  scalar  field  C  to  any  Hilbert 
space  HI,  and  also  generalizing  the  root  of  unity  e^  in  (2)  to  a  unitary  transformation 
U  :  HI  — >■  HI  of  finite  order.  In  particular,  whereas  the  DFT  is  defined  over  £(Zn,  C),  as 
given  in  (1),  the  GST  is  an  operator  over  the  space  of  all  discrete-time,  vector- valued 
TV-periodic  sequences: 

£{Zn,  H)  :=  {X  :  Z  ->  H  |  X[n  +  N]  =  X[n],Vn  £  Z }.  (3) 

Specifically,  the  GST  is  an  operator  from  this  space  into  itself: 
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Definition  1.  Let  U  be  a  unitary  operator  on  El  with  the  property  that  UN  =  I.  The 
Geometric  Sum  Transform  of  any  X  e  ^(Zjv,  El)  is  GST[/(AT)  e  £(7in,M)  whose  nth 
element  is  the  vector: 

N-l 

GSTt/(X)[n]  :=  U~nn' X[n'\.  (4) 

n'= 0 

Note  that  the  order  of  a  unitary  operator  U  is  defined  as  the  least  nonnegative  integer 
K  such  that  UK  =  I.  In  order  to  define  the  GST,  the  order  of  U  must  necessarily 
divide  N. 

Like  the  DFT,  the  GST  (4)  can  be  used  to  estimate  frequency.  While  the 
DFT  determines  the  frequency  of  a  sequence  of  numbers  by  measuring  how  quickly 
they  circle  the  origin  in  the  complex  plane,  the  GST  estimates  the  frequency  of  a 
sequence  of  vectors  with  respect  to  the  action  of  U.  To  be  precise,  if  U  is  a  rotation 
operator  then  GST[/(AT)  can  be  used  to  estimate  the  angular  frequency  of  objects  seen 
rotating  in  the  video  sequence  X.  Meanwhile,  if  U  is  a  cyclic  translation  operator, 
then  GST[/(X)  can  be  used  to  estimate  lateral  velocity. 

To  see  this,  consider  the  toy  example  depicted  in  Figure  2.  There,  N  —  8,  the 
Hilbert  space  is  El  =  L2(M2),  and  the  video  sequence  is  thus  a  member  of  f?(Zg,  L2(M2)). 
The  sequence  of  images  X[n]  is  depicted  in  the  first  row,  where  n  —  (0, . . . ,  7).  It 
turns  out  the  Mercedes-Benz  shape  in  these  images  is  rotating  by  a  factor  of  3(^L) 
radians  in  each  image  frame,  namely  by  135  °  counterclockwise  from  the  previous 
frame.  Suppose  we  want  to  use  the  GST  to  determine  this  fact.  Since  N  —  8,  we  can 
only  apply  the  GST  with  an  operator  U  that  satisfies  Us  =  I.  We  choose  U  to  be 
rotation  by  namely  45°  in  the  counterclockwise  direction. 

To  see  what  the  GST  does  in  this  example,  note  that  the  JJ~nn'  term  in  the 
GST  definition  (4)  fixes  n  and  then  sums  over  different  values  of  n! .  Specifically,  for 
a  conjectured  frequency  n,  the  GST  simulates  taking  a  long  exposure  as  the  camera 
rotates  backwards — in  the  clockwise  direction — at  a  rate  of  n  complete  revolutions 
over  the  N  image  frames.  To  see  this,  the  X[0]  column  of  Figure  2  is  equivalent  to 
holding  the  camera  still,  that  is,  rotation  by  0  radians  per  image.  Meanwhile,  the  A"[l] 
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column  corresponds  to  rotating  the  camera  by  —  ^  radians  per  image.  This  process 
is  continued  for  each  possible  frequency  until  the  camera  has  been  rotated  by  —  7(^) 
radians  in  the  final  column. 

When  the  conjectured  frequency  n  corresponds  to  the  true  rotational  frequency, 
the  two  rates  of  rotation,  both  simulated  and  actual,  cancel  each  other  perfectly. 
This  yields  a  coherent  sum,  resulting  in  a  sharp  image.  For  any  other  frequency,  the 
corresponding  sum  is  at  least  partially  incoherent,  resulting  in  a  blurry  image.  One 
can  thus  determine  the  rate  of  rotation  by  computing  the  sharpness  of  the  sum.  In 
Figure  2,  GST(A"[3])  appears  to  be  the  sharpest  GST  image;  this  confirms  that  the 
Mercedes-Benz  shape  indeed  makes  3  complete  revolutions  over  the  course  of  this 
simulated  video  sequence. 

Meanwhile,  for  the  100-frame  thruster  video  sequence,  15  frames  of  which  are 
depicted  in  Figure  1,  the  resulting  GST  consists  of  100  images.  Rather  than  plot  each 
of  these  images  and  look  for  the  one  which  is  the  sharpest,  we  can  alternatively  plot  a 
sequence  of  their  norms,  as  shown  in  Figure  3.  Intuitively,  the  peak  of  this  norm  plot 
corresponds  to  the  true  rate  of  rotation,  as  we  expect  coherent  sums  to  have  larger 
norms  than  incoherent  ones.  Later  in  this  thesis,  we  formally  show  that  this  is  indeed 
the  case. 

From  the  above  examples,  we  see  that  the  GST  is  useful.  However,  the  question 
now  arises  of  whether  or  not  it  is  practical  to  compute.  Since  n  and  n'  both  have 
N  possible  values  in  (4),  the  computation  of  the  GST  appears  to  require  0(N2) 
matrix-vector  multiplications.  When  N  is  large,  this  can  be  prohibitively  expensive. 
Fortunately,  in  this  thesis  we  show  that  like  the  DFT,  there  exist  a  fast  algorithm  for 
computing  the  GST  that  reduces  this  computational  burden  to  0(N  log  N)  matrix- 
vector  multiplications. 
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1.1  Major  contributions 

The  most  significant  contribution  of  this  thesis  is  the  Fast  Geometric  Sum  Trans¬ 
form  (FGST),  a  fast  algorithm  to  compute  the  GST.  Applied  to  a  video  sequence,  the 
FGST  significantly  decreases  the  amount  of  matrix-vector  multiplications  required  to 
compute  the  GST.  Indeed,  the  FGST  improves  the  GST  computations  required  from 
0(N2)  to  0(N  log  N)  matrix- vector  multiplications.  In  this  way,  we  see  that  the  GST 
mimics  two  important  facets  of  the  DFT:  (1)  it  can  be  used  to  estimate  frequency,  and 
(2)  it  can  be  computed  surprisingly  quickly.  Additionally,  there  are  other  properties 
of  the  GST  that  are  similar  to  those  of  the  DFT,  and  are  proven  in  this  thesis.  For 
example,  the  GST  is  a  linear  operator,  has  an  adjoint,  and  can  formally  be  shown  to 
provide  a  good  frequency  estimate  for  a  certain  class  of  signals.  Further,  we  perform 
a  spectral  analysis  which  shows  that  the  GST  is  invertible  under  certain  conditions, 
and  moreover  shows  that  the  GST  can  indeed  be  regarded  as  a  system  of  DFTs. 

Several  of  these  results  are  novel  contributions  to  mathematical  literature.  In 
particular,  as  far  as  we  know,  our  definition  of  the  GST  (4)  is  the  first  time  any¬ 
one  has  generalized  the  DFT  in  a  way  that  allows  one  to  estimate  rotational  fre¬ 
quency  from  video.  Moreover,  the  proof  of  the  FGST  algorithm  (Theorem  5),  though 
heavily  inspired  by  the  well-known  FFT  algorithms,  is  substantially  more  succinct 
and  straightforward  than  all  other  similar  decimation-in-frequency  arguments  that 
we  could  find  in  the  existing  literature.  Finally,  the  spectral  theory  of  Chapter  V, 
which  formally  shows  that  the  GST  is,  in  fact,  equivalent  to  computing  a  system  of 
DFTs,  is  completely  new. 

1.2  Outline 

Chapter  II  provides  a  summary  of  previous  literature  in  this  area  of  research. 
Since  this  thesis  introduces  the  GST,  we  focus  on  background  theory  about  the  DFT, 
as  these  two  transforms  have  many  of  the  same  properties.  Many  basic  properties 
of  the  GST  are  provided  in  Chapter  III.  Then,  Chapter  IV  includes  a  clear  way  to 
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define  a  generalized  bit  reversal,  and  most  importantly,  provides  the  FGST  algorithm 
and  the  proof  of  its  validity.  A  spectral  analysis  of  the  GST  is  given  in  Chapter  V 
and  we  conclude  in  Chapter  VI  with  more  examples  and  numerical  experimentation 
concerning  the  GST.  The  appendix  includes  Matlab®  code  for  the  implementation 
of  the  FGST,  so  that  this  research  can  be  more  easily  applied  in  future  research. 
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GST(XfOl)  GST(Xfll)  GST(X[2l)  GST(X[3l)  GST(X[4l)  GST(X[5l)  GST(X[6l)  GST(X[7l) 


Figure  2:  This  is  a  GST  rotation  example  of  data  length  N  =  8  applied  to  the 
Mercedes-Benz  shape.  The  first  row  is  the  original  X  sequence.  Each 
column  under  that  initial  row  shows  successively  rotated  images  that  are 
added  together  to  form  the  final  image  in  each  column,  namely  the  GST  of 
a  given  X  [n] . 
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Figure  3:  This  plot  shows  the  £°°  norm  of  the  GST  of  a  thruster  video,  several  image 
frames  of  which  are  depicted  in  Figure  1.  In  this  example,  N  =  100  and  U  is 
rotation  by  Intuitively,  it  makes  sense  that  the  correct  rate  of  rotation 
corresponds  to  the  peak  in  the  plot,  since  at  the  true  rate  of  rotation,  we 
expect  the  image  frames  to  add  coherently,  in  particular  at  their  brightest 
point.  Here,  that  peak  occurs  at  n  =  8,  which  once  converted  to  standard 
units  means  the  plasma  is  rotating  at  a  rate  of  80  Kilohertz. 


II.  Literature  Review 


The  GST  (4)  is  a  new  way  to  generalize  the  DFT  (2).  Indeed,  it  is  a  natural  gen¬ 
eralization  because  of  the  relationship  between  the  DFT  and  the  GST.  Recall  that 
the  FFT  is  a  fast  algorithm  to  compute  the  DFT.  In  this  thesis,  we  provide  the 
FGST,  a  similarly  fast  algorithm  for  computing  the  GST.  In  this  chapter,  we  focus 
on  background  theory  about  the  DFT  alone,  since  the  GST  and  the  DFT  share  many 
of  the  same  properties,  and  no  previous  literature  about  a  GST-like  operator  could 
be  found. 

Three  properties  of  the  DFT  that  are  related  to  the  work  presented  in  this  thesis 
are  (1)  the  DFT  has  an  inverse,  (2)  the  DFT  is  a  linear  transform,  and  (3)  the  FFT 
assumes  periodicity  [11].  In  addition,  the  FFT  computes  the  finite  Fourier  transform 
of  a  length  N  data  set  in  0(N  log  N)  operations  [4],  Contrast  this  difference  in 
efficiency  to  the  0(N'2)  operations  required  to  compute  the  DFT  directly. 

Two  main  reasons  the  FFT  has  a  computing  advantage  over  directly  computing 
the  DFT  is  because  it  applies  known  periodicity  of  the  sine  and  cosine  functions 
and  also  exploits  indices  that  contain  a  factor  that  is  a  power  of  two  [4],  Cooley  and 
Tukey  [3]  showed  how  to  compute  an  FFT  whose  length  is  a  highly  composite  number 
N,  that  is,  an  integer  with  many  prime  factors.  Their  research  showed  that  an  FFT 
of  size  N  =  pip2  . .  .pm.  takes  N(pi  +  p2  +  . . .  +  pm )  operations. 

Often,  the  first  step  in  many  FFT  algorithms  of  length  N  is  to  reverse  the  order 
of  the  bit  representation  of  the  indices  n  —  (0, . . . ,  N  —  1).  Such  bit  reversal  can  be 
found  by  using  Buneman’s  algorithm,  as  presented  by  Walker  in  [14],  For  example, 
the  number  5  has  a  bit  representation  of  0101,  which  when  reversed  becomes  1010 
and  corresponds  to  the  number: 

/3( 5)  =  1(23)  +  0(22)  +  1(21)  +  0(2°)  =  10. 

This  description  of  bit  expansion  and  this  example  are  described  in  more  detail  in  [1], 
In  Chapter  IV,  we  generalize  this  notion  of  bit  reversal  to  arbitrary  n-ary  expansion. 
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Although  all  FFT  algorithms  incorporate  bit  reversal,  each  algorithm  can  be 
computed  using  either  decimation  in  time  or  decimation  in  frequency.  The  1965 
Cooley- Tukey  algorithm  is  an  example  of  decimation  in  time,  that  is,  where  each 
computation  group  gets  bigger  than  the  previous  one  [11].  On  the  other  hand,  the 
1966  Gentleman-Sande  algorithm  [5]  is  an  example  of  decimation  in  frequency,  that  is, 
where  each  computation  group  is  smaller  than  the  previous  level.  These  two  different 
computation  methods  are  best  visualized  with  well-known  FFT  “butterfly  diagrams,” 
as  shown  in  [1],  [2],  and  [14], 

According  to  Heideman  [6],  in  1805  preceding  Fourier’s  introduction  of  sinu¬ 
soidal  expansion,  Gauss  created  an  algorithm  that  closely  resembles  the  FFT,  but  his 
algorithm  did  not  quantify  computational  complexity  to  the  N  log  N  operations  FFT 
algorithm  we  now  use.  Gauss’s  algorithm  computed  the  coefficients  of  a  finite  Fourier 
series,  and  was  similar  to  a  decimation  in  frequency  algorithm  that  has  been  modified 
for  a  real  data  sequence  [6].  Note,  the  FGST  algorithm  we  present  in  Chapter  IV 
is  also  an  example  of  decimation  in  frequency,  where  the  initial  input  length  value 
is  factored  and  then  computation  groups  are  continually  broken  down  into  smaller 
groups. 

Now,  more  than  two  hundred  years  after  Gauss’s  work,  there  have  been  many 
different  FFT  algorithms  created.  In  [13],  Van  Loan  presents  a  summary  of  many 
popular  FFT  algorithms,  and  categorizes  them  based  on  the  size  of  each  framework. 
That  is,  the  frameworks  are  distinguished  by  how  an  FFT  of  length  N  is  broken  into 
computation  groups.  For  instance,  there  is  a  class  of  algorithms  specifically  created 
for  computing  the  FFT  of  a  length  N  =  2k  signal  for  k  E  N.  Several  examples  of 
these  algorithms  are  the  Cooley- Tukey  algorithm,  the  Stockham  algorithm,  and  the 
Pease  algorithm,  all  presented  in  [13]. 

Variations  of  each  algorithm  are  also  presented  in  [13] ,  including  ways  of  splitting 
the  computational  groups,  such  as  Radix-2,  Split-Radix,  and  Mixed-Radix.  A  Radix- 
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2  framework  is  applied  to  a  signal  of  length  N  =  2k,  but  a  Mixed- Radix  framework 
is  applied  to  a  signal  of  length  N  —  pi  ■  ■  ■  pk  where  each  p*  is  a  factor  of  N. 

Most  of  the  FFT  algorithms  we  have  mentioned  up  to  this  point  are  similar 
to  the  FGST  algorithm  we  introduce  in  Chapter  IV.  However,  there  are  some  FFT 
algorithms  that  exactly  follow  the  approach  we  pursue  in  this  thesis.  For  example,  in 
1968,  Rader  created  a  way  to  compute  an  FFT  of  prime  length  by  treating  the  signal 
like  a  cyclic  convolution  [10].  Then,  in  1999,  Mohlenkamp  found  a  fast  transform  for 
spherical  harmonics,  that  is,  problems  where  the  Fourier  exponential  functions  are 
generalized  to  functions  over  the  sphere  S 2  in  M3  [9].  More  recently  in  2006,  Rokhlin 
and  Tygert  presented  a  fast  algorithm  for  forward  and  inverse  spherical  harmonic 
transforms  [12]. 
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III.  Basic  properties  of  the  GST 

In  this  chapter  we  introduce  and  prove  several  of  the  basic  properties  of  the  GST  (4). 


3.1  £(Zn,HI)  is  a  Hilbert  space 

Recall  the  GST  (4)  is  defined  over  the  space  £( Zjv,HI)  defined  in  (3).  We  now 
verify  that  this  space  is  indeed  a  Hilbert  space  with  inner  product: 

N-l 

(X,Y)«z„m  !•>])„.  (5) 

n=0 

To  verify  this,  we  just  need  to  show  that  (. X ,  Y)^zN,m)  satisfies  the  four  properties 
of  an  inner  product  space.  First,  for  any  X,Y,  Z  e  £(1in,  HI): 

N-l 

(X  +  Y,  Z)e{ZNM)  =  ^<X[n]  +  Y[n],  Z[n])u 

n= 0 
N-l 

=  ^((X[n],Z[n]U  +  (Y[n],Z[n])K) 

n= 0 

N-l  N-l 

=  ^2(X[n\,  Z[n]) H  +  ^2(Y[n],  Z[n])m 

n= 0  rc=0 

=  (X,  Z)i(zN, H)  +  (Y,  Z)^zN, H). 


Second,  for  any  a  G  C  and  any  X,  Y  e  £(Zat,  HI): 

N-l  N-l 

(XiOtY)^iNjg)  =  ^{X[n],«y[n])i  =  y^(X[n],  aY[n])  H  =  a(X,  Y)g(zN,M)- 

n= 0  n— 0 
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Third,  for  any  X,Y  e  £(ZN,M)\ 


N-l 

(X,Y)e{zN,  h)  =  ^<A>],1>])h 

n= 0 
N-l 

=  E<rw,^w>i 

n= 0 
AT— 1 

EO'KAW. 

n=0 


Lastly,  for  A"  e  £(Zjv,EI)  with  A”  7^  0,  we  have  Ai[n]  7^  0  for  some  n  and  so: 

N-l  N-l 

(X,X}t^,H)  =  ^(A[nU[n])H  =  V  ||.Y[n]||ji  >  0. 

n=0  n= 0 

Therefore  (5)  is  indeed  an  inner  product  space  on  (3),  as  claimed. 


3.2  GSTu  is  a  linear  operator 

To  better  understand  the  GST  operator,  we  first  show  that  GST;y  is,  for  any 
unitary  operator  U,  a  well-defined  linear  operator  from  ^(Zjv,HI)  into  itself.  That  is, 
the  domain  of  the  definition  of  GSTy  is  £(Zjv,EI).  Then  for  any  X  e  ^(Zjv,  H),  the 
fact  that  GST[/(A")  is  TV-periodic  follows  from  the  assumption  that  UN  =  I: 

N-l 

GST u{X)[n  +  N}  =  ^  U~n'{n+N)X[n'} 

n'= 0 

N-l 

=  J2U~nn'  {UN)~n'  X[n'} 

n'= 0 

JV-1 

=  J2  U-nn'X[n'} 

n'= 0 

=  GSTt/(X)[n], 
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Having  this  fact,  we  next  show  that  GSTV  is  linear.  Let  X]}  X2  be  vectors  of  length 
N  such  that  both  X1}X2  E  ^(Ztv,  HI).  Then  for  all  n  E  Z^v  and  scalars  au,  a2  E  C,  we 
fix  n  and  have: 

N—l 

GSTv(aX  +  pY)  [n]  =  Y  U~n'n(aX[n ']  +  pY[n ']) 

n'= 0 

N—l 

=  J2  (aU~n'nX[n']  +  pU~n'nY[n']) 

n'= 0 

N—l  N—l 

=  a  U-n'nX[n'}  +PY  U~n'nY[n'} 

n'= 0  n'= 0 

=  aGSTu(X)[n]  +  /3GSTc/(F)[n] 

=  aGST V(X)  +  /5GSTf/(r)l  [n\. 

Since  this  is  true  for  all  n  E  ZN,  a,  /3  E  C  and  X,  Y  E  £(Z N ,  H),  then  GST  is  a  linear 
operator. 

As  discussed  in  the  introduction,  the  GST  is  a  generalization  of  the  DFT.  More¬ 
over,  one  of  the  well-known  properties  of  the  DFT  is  that  conjugating  a  translation 
operator  by  the  DFT  results  in  a  modulation  operator.  Specifically,  defining  the 
operator  Tfe  :  ^(Zjv,HI)  — >  ^(Zjv,HI),  (Tfcx)[n]  :=  x[n  —  k ],  it  is  easy  to  show  that 
DFT(Tfca;)[n]  =  e“^DFT(a:)[n]  for  any  k,n  E  Z.  If  we  generalize  this  notion  of 
translation  to  the  operator  Tfc  :  £(ZN,M)  — y  £(ZNlM),  (TfcA")[n]  :=  X[n  —  k]  then, 
making  the  substitution  n"  —  n'  —  k  gives  a  similar  result: 

JV-1 

GSTlf(TfcX)[n]  =  Y  U~nn'X[n'  -  k\ 

n'= 0 

N—l 

=  U~kn  Y  U~nn"x[n "] 

n"= 0 

=  U~knGSTu  (X )  [n] . 
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3.3  The  adjoint  of  the  GST 

The  GST  also  generalizes  the  well-known  characteristics  of  the  adjoint  of  the 
DFT,  namely  the  operator  DFT*  :  £(Z jv)  — y  I^Ln)  where: 

N- 1 

_  .  .  .  r  X  >  27rinn  r  1-1  ,  , 

DFT  (x)[n\  =  ^e  N  x[n].  (6) 

n'=0 

In  particular,  the  only  difference  between  DFT  (2)  and  its  adjoint  (6)  is  the  sign  of 
the  exponential.  We  now  show  a  similar  result  holds  for  the  GST. 

Theorem  2.  GST^  =  GST^*. 

Proof.  It  suffices  to  show  that  for  any  two  vectors  X ,  Y  e  £(Zat,  H): 

(X,  GST*v(Y))nZNM)  =  (X,  GST u*(Y))£{ZnM). 

By  definition  of  the  adjoint  and  the  inner  product  (5)  on  £(Zjv,EI),  we  have: 

(X,GST^(F)),(Zjv>h)  =  (GST^(X),F)t(Zjv,H) 

N- 1 

=  ^(GST„(.Y)[n].F[n]) 

\  /  M 

n= 0 

JV— 1  JV-1 

n=0  n'= 0 
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Interchanging  sums  and  using  the  fact  that  U  is  unitary  gives: 


N—l  N- 1 

(X,GSTu(Y))e{ZNtM)  =  EE(  U~nn' X  [n'\ ,  Y [n]  ^ 

n= 0  n/=0 
N—l  N—l 

=  EE(xM-(t7""”')*yw)H 

n'= 0  n= 0 

V-l  V-l 

=  E(*M>Et,””'3>])II 

n'= 0  n=0 

JV-1 

=  E(-’fM.GST„-,(y)K]) 

n'= 0 

=  (x,GSTt7.(y)\ 

\  /  ^(Zjv  ,H) 


Despite  these  many  similarities  between  the  DFT  and  the  GST,  they  are  not 
completely  identical.  Indeed,  the  DFT  is  invertible  with  DFT-1  =  TdFT*.  However, 
this  is  not  always  the  case  for  the  GST.  For  example,  when  U  —  I,  GST/  produces 
N  copies  of  the  sum  of  the  entries  of  X ,  and  so  GST /  is  clearly  not  invertible.  It  will 
be  shown  with  more  in-depth  spectral  analysis  in  Chapter  V,  GST//  is  only  invertible 
for  certain  choices  of  U . 


3-4  Estimating  frequency  with  the  GST 

Applications  discussed  in  Chapter  VI  have  HI  =  L2(M2),  where  the  most  impor¬ 
tant  property  of  the  GST  is  its  ability  to  estimate  rate  of  rotation.  To  be  precise, 
the  dominant  frequency  of  a  complex-valued  signal  x  is  usually  taken  as  the  location 
of  the  highest  peak  in  its  DFT.  The  GST  has  a  similar  property:  given  video  of  an 
object  rotating  at  a  uniform  rate,  the  next  result  shows  how  this  rate  of  rotation  can 
be  computed  by  finding  the  “peak”  of  the  GST. 

Theorem  3.  If  X  G  ^(Zjv,HI)  has  the  property  that  there  exist  no  G  Z  such  that 
X[n+ 1]  =  Un°X[n ]  for  alln  G  Z,  thenno  =  argmaxn  ||GST//(V)[n] ||,  where  ||  -||  may 
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be  any  U -invariant  norm,  that  is,  any  norm  that  has  the  property  that  \\Ux\\  =  ||a: 
for  all  x  G  HI. 


Proof.  The  given  assumption  is  equivalent  to  having  X[n]  =  Unn°X[0 ]  for  all  n  e  Z. 
Using  this  property,  we  have: 


GSTc/(X)[n]|| 


N- 1 

Y  U~nn'X[n'} 

n'= 0 


N- 1 

Y  U-nn'Unon'X[ 0] 

n'= 0 


N- 1 

=  Yun'{no~n)x[°]  ■ 

n'= 0 

By  the  triangle  inequality  and  the  U-invariance  of  the  norm: 


N-l 


Y  Un'{n°-n)X{ 0] 


n'= 0 


N-l 

n'= 0 


Un'(n  0-»)y[0] 


N-l 


^||A'[0]||=JV||A'[0]||, 


n'= 0 


Therefore  ||GST[/(A")||  has  an  upper  bound  of  iV||A"[0]||.  Moreover,  this  upper  bound 
is  achieved  at  n  =  n0: 


||GSTc/(X)[n0]||  = 


N-l 


Y  Un'{n°-no)X[ 0] 


n'= 0 


N-l 


E  A'i°i 


n'= 0 


=  A||A[0]||. 


Hence,  iV|| AU[0]  ||  is  the  maximum  value  achieved  by  any  choice  of  n,  and  so  no  = 
argmaxn  ||GSTti(A^)[n]||,  as  claimed.  □ 

Note  that  no  is  not  necessarily  unique.  For  example,  when  U  =  I  the  GST  is 
constant  and  so  any  choice  of  n0  is  equal  to  arg  maxn  ||GSTC/(A")[n]||.  In  summary, 
the  previous  result  shows  how  a  GST,  once  computed,  can  provide  useful  information 
about  the  rate  of  rotation  (or  translation)  in  a  given  sequence  of  vectors.  However, 
as  we  shall  see  in  the  following  chapter,  the  GST  can  be  expensive  to  compute,  both 
in  time  and  memory,  whenever  N  or  the  dimension  of  HI  is  large.  This  raises  the 
question:  is  there  a  fast,  memory  efficient  algorithm  for  computing  a  GST?  In  the 
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next  chapter,  we  show  that  such  an  algorithm  indeed  exists.  In  particular,  since  the 
GST  is  a  generalization  of  the  DFT,  we  look  to  generalize  the  well-known  decimation- 
in-frequency  fast  Fourier  transform  algorithm  to  the  GST  setting. 
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IV.  A  fast  algorithm  for  computing  the  GST 

In  the  previous  chapter,  we  discussed  how  the  GST  is  a  useful  tool  for  estimating 
frequency.  Note  however,  that  computing  the  GST  directly  requires  N'2  matrix-vector 
multiplications,  that  is  N  multiplications  for  each  of  the  N  distinct  choices  of  n  in 
(4).  In  this  chapter,  we  will  present  a  fast  algorithm  for  computing  the  GST  which 
generalizes  the  way  in  which  the  FFT  quickly  computes  the  DFT. 

As  mentioned  in  the  introduction,  Cooley  and  Tukey’s  FFT  algorithm  takes 
advantage  of  the  property  that  an  DFT  of  a  given  nonprime  size  can  be  successively 
broken  into  smaller  DFTs  of  prime  order.  When  N  is  not  prime,  the  GST  behaves 
similarly.  Indeed,  if  N  is  divisible  by  P  G  N,  then  an  Appoint  GST  can  be  computed 
in  terms  of  P  GSTs  of  length  To  see  this,  note  that  making  the  substitution 
vl  =  Pq  +  p  gives: 

N—l 

GSTc/(X)[n]  =  J2  U-nn'X[ri] 

n'= 0 

P-1  f-1 

U~(Pq+p)nX[Pq  +  p] 

p= o  <j=0 

p-i  f-i 

=Eu~p"E(i ’prnx[p<i+vi 

p= 0  q= 0 

By  the  definition  of  GST,  we  have: 

p-i 

GSTf/(X)[n]  =  Y,U~pnGSTuP(X[P(0  :  f  -  1)  +  p])H 

p= o 

where  the  notation  A"[P(0  :  ^  —  1)  +  p]  denotes  the  subsequence  of  X  whose  indices 
lie  in  the  coset  {Pq  +  p  :  q  =  0, . . . ,  ^  —  1}. 

For  example,  a  16-point  GST  can  be  broken  into  two  GSTs  of  size  4,  or  alterna¬ 
tively  eight  GSTs  of  size  2,  or  even  two  GSTs  of  size  8.  This  property  also  holds  for 
the  FFT.  In  fact,  in  the  literature,  a  size  2k  FFT  for  some  fc  £  N  is  referred  to  as  a 
radix- 2  FFT.  Meanwhile,  a  FFT  of  size  N  where  N  can  be  factored  into  non-primes  is 
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sometimes  referred  to  as  a  mixed-radix  FFT.  Later  in  this  chapter  we  show  that  the 
prime  factorization  is  the  particular  factorization  which  makes  our  FGST  algorithm 
perform  at  its  best. 

Moreover,  in  such  algorithms,  it  is  not  only  the  factors  themselves  that  matter, 
their  order  is  also  significant.  Indeed,  a  fundamental  part  of  many  FFT  computations, 
namely  the  bit  reversal  operation  defined  below,  is  defined  in  terms  of  this  ordered  list 
of  factors.  Note  that  we  are  using  the  term  “bit  reversal”  to  imply  an  n-ary  represen¬ 
tation  reversal.  For  example,  when  N  =  2k  for  some  k  €  Z,  the  binary  representation 
of  any  number  in  the  range  (0, . . . ,  N  —  1)  is  length  k  with  each  “significant  digit” 
having  value  either  0  or  1.  A  bit  reversal  reverses  the  order  of  these  digits.  More 
generally,  for  any  N,  the  number  of  prime  factors  of  N  determines  the  length  of  the 
n-ary  representation,  where  each  significant  digit  is  determined  by  the  prime  factors 
themselves.  To  be  precise,  we  have  the  following  definition. 

Definition  4.  Let  N  €  N  and  let  V  =  {P0  •  •  -Pj_i}  be  an  ordered  factorization  of 
N,  that  is  N  —  n/=o  where  each  Pj  is  a  positive  integer  greater  than  1.  We  define 
the  partial  product  sequence  as  J\f  =  {N0, . . . ,  Nj},  where  for  any  j  —  (0, . . . ,  J), 

j- i 

JVf  •  II  Pf 

j'=j 

which  gives  Nj  =  1,  and  we  also  observe: 

Nj  —  PjNj+1,  Vj  =  (0, . . . ,  J  —  1).  (7) 

For  each  j  —  (1, . . . ,  J)  we  define  the  bit  reversal  f3j  as  a  permutation  on  {0, . . . ,  — 

1}.  This  definition  is  recursive.  Specifically,  we  let  /?i(n)  :=  n  for  all  n  =  (0  . . . ,  P0  —  1). 
Meanwhile,  for  any  j  =  (1 . . . ,  J  —  1)  we  let: 

Pj+i (pjn  +  P)  ■=  -w~P  +  Pj(n),  Vn  =  (0, . .  -  1),  p  =  (0, . . . ,  Pj  —  1).  (8) 
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For  example,  suppose  N  =  30  and  let  P0  =  2 ,  Pi  —  3,  and  P2  =  5.  The 
partial  products  for  this  ordered  factorization  of  N  are  N0  =  30,  Ni  =  15,  N2  =  5, 
and  N3  =  1.  So,  using  (8),  n  =  27  has  the  n-ary  expansion  27  =  2(1)  +  2(5)  + 
1(15)  =  2IV3  +  2N2  +  1  Ni.  Then,  to  compute  j33( 27),  we  reverse  the  bits  such  that 
/33( 27)  =  l-^  +  2-^-  +  2-^  =  1(1)  +  2(2)  +  2(6)  =  17.  However,  computing  the  bit 
reversal  in  this  way  is  cumbersome.  An  equivalent  but  faster  method  for  computing 
the  bit  reversal  is  to  use  Kronecker  tensor  products: 

N 

Pj+i  ■=  Pi  ®  Ip,  +  Wl jv  ®  (0  :  P3  -  1).  (9) 

The  <8)  symbols  above  represent  the  Kronecker  product  operations,  also  known  as 
matrix  direct  products,  and  the  the  1^  denotes  an  N  x  1  vector  of  ones.  This  method 
of  computing  bit  reversal  is  also  implemented  in  MATLAB  code  given  in  the  appendix. 
Returning  to  the  N  =  30  example,  we  will  now  show  how  to  recursively  compute  bit 
reversal  using  (9): 

0 
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The  binary  bit  reversal  discussed  in  the  introduction  exploited  the  characteris¬ 
tics  of  binary  bit  flipping.  The  n-ary  bit  reversal  given  in  Definition  4  generalizes  this 
operation  to  any  N.  Further  note  that  these  bit  reversals  are  well-defined  permuta¬ 
tions.  We  can  see  this  by  induction  on  j.  By  definition,  Pi(n)  =  n,  and,  thus,  j3\  is 
one-to-one.  Now  suppose  (3j  is  one-to-one.  To  see  that  (3j+i  is  also  one-to-one,  assume 
for  some  ki,  k2  =  (0, . . . ,  —  1)  that  j3j+i(k\)  =  ftj+\(k2).  Write  k\  =  PjUi  +  p\ 

and  k2  =  Pj‘n2  +  p2  for  some  rii,n2  =  (0, . . . ,  ^  —  1)  and  pi,p2  =  (0, . . . ,  Pj  —  1),  we 
have: 


§:P\  +  PiM  =  Pj+iiPjTi!  +pi)  =  Pj+1(Pjn2  +  p2)  =  j$-p2  +  Pj(n2). 

Subtracting  gives  -p2)  =  Pj{n2)  -  Pj(ni).  Since  pj(ni),  Pj(n2)  G  [0, . . . ,  ^  -  1] 
this  implies  jj-(p\  —  p2)  =  0.  Thus,  p\  =  p2,  and  so  f3j(ri\)  =  /3j(n2),  implying  by 
the  inductive  hypothesis  that  ri\  =  n2.  Thus,  k\  =  k2  and  so  /3j+ 1  is  one-to-one,  as 
claimed. 

In  the  next  result,  we  show  how  this  notion  of  bit  reversal  can  be  used  in  a  fast, 
memory  efficient  algorithm  for  computing  the  GST.  In  short,  this  result  shows  that 
the  decimation-in-frequency  method  for  computing  an  FFT  generalizes  to  the  GST 
setting. 
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Theorem  5.  The  GST  of  any  X  e  £( ZN,W)  can  be  computed  by  the  following  Fast 
Geometric  Sum  Transform  (FGST )algorithm: 


(00)  let  N  e  N,  JeN  with  1  <  J  <  N 

(01 )  let  [n]  =  X[n]  for  all  n  =  (0, . . . ,  N  -  1) 

(02)  for  j  —  (0, . . . ,  J  —  1) 

( 03)  for  k  =  (0, . . . ,  Nj+i  -  1) 

(04)  forl  =  (  0, 

(05)  forp  =  (0,. . .  ,Pj  -  1) 

(06)  X^[NJ+1(Pjl+p)  +  k\ 

=  E(^ (^>[JVj+iWI+4)+fc] 

g=o 

(07)  end 

(08)  end 

(09)  end 
(10)  end 

Then  for  all  n  —  (0, . . . ,  N  —  1)  ; 

X^J\n\  =  GSTj/(X)[/3j(n)].  (10) 


Proof.  For  brevity  of  notation,  let  IV  e  N  and  J  G  N  with  J  <  N,  then:  GST[/(X)  := 
GST(X,  U).  For  any  j  =  (0, . . . ,  J)  we  claim  it  suffices  to  show  the  values  of  the  jth 
iteration  of  the  GST  algorithm  are  given  by: 

X^[Nja  +  b}  =  GST(X[Nj(0  :  ^  -  1  )  +  b],  UNi)[^(a)\  (11) 

for  all  a  =  (0, . . . ,  —  1)  and  b  =  (0, . . . ,  Nj  —  1),  where  /3j  is  the  -^-ary  bit  reversal 

permutation  on  {0, . . . ,  jf-  —  1}  given  in  Dehnition  4.  Indeed,  if  (11)  holds  for  all 
j  =  (0, . . . ,  J),  then,  in  particular,  when  j  =  J  we  have  a  —  (0, . . . ,  N  —  1)  and  b  =  0. 
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Letting  a  —  n,  (11)  becomes: 


X^[n]  =  GST(X, U)[(3j(n)}. 

We  prove  (11)  by  induction  on  j.  When  j  =  0,  we  have  a  =  0  and  b  =  (0, . . . ,  N  —  1) 
where  the  right  hand  side  of  (11)  becomes: 

GST(A'[f>],/)][/3„(0)]  =  A[6]  =  A<°>[6]. 

Hence  (11)  holds  for  j  =  0.  Now,  assume  (11)  holds  for  a  given  j  e  {0, . . . ,  J}.  We 
want  to  show  (11)  holds  for  j  +  1,  that  is: 

a :{i+1)[Nj+1(pji  +  P)  +  k] 

=  GST(X[JV,+1(0  :  ^  -  1)  +  k\,UN‘»)[pj+l(Pjl  +  V)]  (12) 

where  l  =  (0  . . . ,  ^  —  1),  p  —  (0, . . , ,  Pj  —  1),  and  k  —  (0, . . . ,  Nj+1  —  1).  The  left  hand 
side  of  (12)  is  given  by  Line  (06)  of  the  algorithm.  To  see  that  this  equals  equals  the 
right  hand  side  we  expand  it: 


GST(X[1V,+1(0  :  ^  -  1)  +  k],UNj+1)\j3j+1(Pjl  +  p)] 

-Ji _ i 

Ni+ 1 

=  (UNj+1)-cPj+l('Pjl+p^ X[Nj+ic  +  k\.  (13) 
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Letting  c  =  Pja  +  q  where  a  =  (0,  —  1)  and  q  —  (0, . . . ,  Pj  —  1),  (13)  becomes: 


GST(X[7V.+1(0  :*£--!)  +  fc],  +  p)} 


vj+i 


-i 


EE  (£/^-+i)-(p^«)^+i(pi*+P)X[iVi+1(Pja  +  q)  +  k] 

q= 0  a=0 

Pj- 1  ^_1 

^  ([/^+i)-^+i(p^+p)  ^  (^)-a^+1(p'z+^X[lVia  +  (Nj+1q  +  k)] 

q= 0  a=0 

Pj~ 1 

(^/Arj+i)-9/3j+i(pj;+p) 

g=0 

■  GST(A'[JV,(0  :  f  -  1)  +  ( Nj+lq  +  k)],UK‘)\fij+i(Pjl  +  p)]. 


(14) 


By  definition,  (3J+\  (P3l  +  p)  —  -^-p  +  Pj(l)-  It  is  also  simple  to  show  that  the  GSTs  in 
(14)  are  ^-periodic.  Thus: 


GST(A[A1,(0  :  f  -  1)  +  (JV,-+l9  +  k)],  f/^JIft+i  W  +  P)\ 

=  GST(A[JV,(0  :  f  -  1)  +  (N1+1q  +  k)}, PK')[f  P  +  ft  W] 

=  GST(A'[Af,(0  :  f  -  1)  +  (JVJ+l9  +  *)],  UN’)\J3j(l)] 

=  X<J>  [Njt  +  NJ+1q  +  k]  (15) 


where  the  final  equality  follows  by  the  inductive  hypothesis  (11).  Then  substituting 
(15)  into  (14)  gives: 


GST(X[Nj+1(0  :  ^  -  1)  +  k\,UNj+1)[Pj+1(Pjl  +  p)\ 

Pj- 1 

=  (UNj+1)-qpj+^Pjl+p)X^[Njl  +  Nj+lq  +  k] 

5=0 

Pj~  1 

=  £  ([/^+i)-^+i(^HP)xb)[Arj+1(p./  +  g)  +  fc]  (16) 

5=0 

which  is  Line  (06)  of  the  algorithm.  Thus,  (11)  holds  for  all  j  —  (0, □ 
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The  computational  requirement  to  run  a  GST  of  size  Pj  is  only  Pj  places  in 
memory.  Matlab  code  to  compute  the  FGST  provided  in  the  appendix  advanta¬ 
geously  uses  in-place  computation  for  memory  efficiency.  Note  that  at  the  heart  of 
the  algorithm  of  Theorem  5  lies  a  GST  of  size  Pj.  Indeed,  the  quantity  computed  in 
Line  (06)  can  be  rewritten  as: 

Pj- 1 

{UNj+1 ) -9di+1  (pjl+p)x(i)  [jVj+1  (Pjl  +  q)  +  k \ 

q= o 

P,-1 

=  (UNi+1)~q{^P+MX{j)[Nj+1(Pjl  +  q)  +  k] 

q= o 

Pj_1  - 

=  ( UPj  )~qpU~Nj+iq^^ X^[(Njl  +  k)  +  Nj+1q}.  (17) 

q= o 

which  is  precisely  the  pth  value  of  GST  of: 

{U-Nj+iqPj{i)xti)  [(Njl  +  k)  +  Nj+lq]}* o1- 

That  is,  the  algorithm  of  Theorem  5  indeed  breaks  a  single  large  GST  computation 
into  many  smaller  ones. 

To  see  that  the  FGST  is  indeed  faster  than  a  direct  computation  of  (4),  note 
that  for  any  fixed  j  =  (0, . . . ,  J  —  1),  Lines  (03)-(09)  involve  (Nj+i)(jk)(Pj)(Pj)  = 
NPj  matrix- vector  multiplications.  Thus,  the  total  matrix- vector  multiplications  is 
N  P]  ■  I11  the  case  when  N  =  PJ ,  that  total  becomes  PNlogp  N  matrix- vector 
multiplications.  For  a  nonprime  N,  this  is  less  than  the  N2  such  multiplications  used 
in  a  direct  computation  of  (4). 

Note  that  the  FGST  algorithm  of  Theorem  5  depends  on  a  given  ordered  fac¬ 
torization  of  N,  as  seen  in  Definition  4.  From  the  point-of-view  of  minimizing  compu¬ 
tation,  the  optimal  factorization  of  N  is  its  prime  factorization.  Indeed,  if  the  factors 
of  N  are  taken  to  be  {Pj}jZo  where  Pj0  =  QoQi  is  not  prime,  then  a  more  factored 
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Figure  4:  A  comparison  of  the  GST  computational  time  using  the  direct  method 
versus  using  the  FGST,  for  length  N  ranging  from  1  to  64. 

algorithm  is  faster;  we  have: 

j-i  J- l 

Nj2PJ>N(Qo  +  Qi  +  J2Pi) 

3=0  3=0 

i¥=io 

since  Q0  +  Qi  <  Q0Qi  =  Pjo. 

Matlab  code  implementing  the  FGST  algorithm  with  a  full  prime  factorization 
is  given  in  the  appendix.  When  N  itself  is  a  prime  number,  the  direct  method  of 
computing  the  GST  takes  the  same  amount  of  time  as  the  FGST,  both  using  N2 
matrix- vector  multiplications.  In  fact,  from  Figure  4  it  is  clear  that  the  FGST  is 
always  at  least  as  fast  as  the  direct  method,  and  is  often  significantly  faster. 

Note  that  here  we  are  measuring  computational  complexity  solely  in  terms  of 
the  number  of  matrix-vector  multiplications.  The  true  computational  cost  depends 
on  the  expense  of  each  such  multiplication,  which,  in  turn,  depends  on  the  nature  of 
U.  We  discuss  these  ideas  in  greater  detail  in  the  following  chapter. 


V.  Spectral  analysis  of  the  GST 

In  the  previous  chapter,  we  gave  a  fast  algorithm  for  computing  the  GST  that  was 
inspired  by  the  classical  FFT  decimation-in-frequency  approach.  In  this  chapter, 
we  show  this  similarity  is  not  a  coincidence;  the  GST  can  be  regarded  as  a  system 
of  DFTs,  provided  we  know  the  eigenvectors  and  eigenvalues  of  U.  To  be  precise, 
assume  El  is  M-dimensional  and  note  that  since  U  is  unitary,  then  it  is  normal.  Then, 
by  Shur’s  Theorem  U  can  be  unitarily  diagonalized  as  U  =  VDV*,  where  D  is  a 
diagonal  matrix  and  V  is  unitary. 

Since  U  is  unitary,  we  further  know  the  diagonal  entries  of  D,  namely  the 
eigenvalues  of  U ,  are  unimodular,  that  is,  have  modulus  one.  Moreover,  since  UN  = 
/,  we  know  these  eigenvalues  are  actually  Nth  roots  of  unity.  Specifically,  letting 
{Kn}m=o  denote  the  eigenvalues  of  U,  we  have  A^  =  1  for  all  m.  Letting  vrn  denote 
the  eigenvector  corresponding  to  Am,  and  writing  Xm  =  e  Nm }  we  have  the  following 
result. 

Theorem  6.  For  any  eigenvector  vm  of  U  with  corresponding  eigenvalue  e— ^ , 

N—l 

.  ,  s  r-  X — ■>  ‘2'K\Q rn.im1  ,  _  . 

{vm,GSTu(X)[n])  =  J^e - -  {vm,X[ri\). 

n/=0 

Proof.  By  the  definition  of  the  GST  given  in  (4), 

,  N—l  v  N- 1 

(vm,GSTu(X)[n])  =  U~nn' x[n'}\  =  ^(um,  U~nn'X[ri]). 

'  n'= 0  n'= 0 

Then  since  U*  =  t/-1,  this  becomes: 

N—l  N—l 

J2(vm,U-nn'X[n'})  =  J2(unn'vm,X[n’}) 

n'= 0  n'= 0 

N—l 

E,  27ri  emnn' 

(e  N  vm,X[n']) 

n'= o 
N—l 

E_  27ri em  nn' 

e  K  (vm,  X  [n  ]).  □ 

n'= 0 
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We  now  use  Theorem  6  to  investigate  the  invertibility  of  the  GST.  In  particular, 
applying  this  result  to  U _1  gives: 


N-l 


(vmi  GST£f(y)[n])  =  (u^GST^y)^])  =  ^  Y[n"}). 


n"= o 


Let  Y  =  GST;y(A")  such  that: 


N—l 


(vm,  GST^(GSTf/(A"))[n])  =  ]T  e^^^ST^AJK]) 

n"= 0 

N-l  n  N-l 

E2-7riflmnn  V — -y  2-K\Qrnn  n  . 

e  N  e  "  (vm,(X)[n'}) 

n"= 0  n'= 0 

N—l  N-l 

=  J2(Vm,(X)in  1)  e' 


2iri6m  (n—n)nn 
N 


(18) 


n'= 0 


n"= 0 


By  the  Geometric  Sum  Formula, 


N-l 

A 

n"= 0 


27ri077i  (n  —  n,)n,> 

e  ^ 


AT, 


2ni0m  (n  — n/) 

e  jv  =  1 


N, 


9m(n-n') 


| _ e2ni6rn(n—n/)  2-jTi6Tn(n—nr) 

2nidrn  (n  —  n/)  1  ^  /" 

1— e  N 


N 


e  z 


o,  ‘My"')  i  z 


(19) 


Then  substituting  (19)  into  (18)  gives: 


iV-l 


(nm,GST^(GSTlf(X))[n])  =  IV  A  A...  ax 


n'=0 

N\9m(n- n') 


N-l 


N(vm, 


n'= 0 

N\9m(n-n') 


Thus,  we  see  that  the  mth  component  of  A:GST^(GSTC/(A"))  is  the  like  compo¬ 
nent  of  the  sum  of  the  entries  of  X  whose  indices  lie  in  a  certain  coset  of  Z^v-  When 
each  0m  is  relatively  prime  to  N,  each  of  these  cosets  only  consists  of  a  single  index, 
meaning  the  GST  is  invertible: 
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Corollary  7.  If  9m  is  relatively  prime  to  N  for  all  m,  then  the  GST  is  invertible  with 

nrr-l  _  1  cinryi* 

JJ. 

If,  on  the  other  hand,  even  a  single  9m  shares  a  common  factor  with  N,  then 
this  summing  process  destroys  some  of  the  frequency  information  in  the  sequence 
{(vm,X[n])}%-0\  This  implies  GST ^GSTu,  and  so  GSTf/,  is  not  invertible  in  such 
cases. 

The  DFT  can  be  viewed  as  a  special  case  of  the  GST  that  is  invertible.  To  be 
precise,  let  H  =  C,  that  is,  N  —  1,  and  let  U  be  the  act  of  multiplying  by  e^,  an 
operation  whose  only  eigenvalue  is  e^,  that  is,  6\  =  1.  Since  9i  is  relatively  prime  to 
N,  Corollary  7  states  that  the  DFT  is  invertible  with  DFT-1  =  TjqjrT*. 

However,  for  the  angular-frequency-determination  problem  that  motivated  this 
work,  the  relatively  prime  condition  of  Corollary  7  does  not  hold,  meaning  the  corre¬ 
sponding  GST  is  not  invertible.  Indeed,  letting  U  be  a  rotation  operator,  the  constant 
vector  v0  is  an  eigenvector  for  U  with  eigenvalue  equal  to  1,  and  so  90  =  0,  which  is 
not  relatively  prime  to  N.  In  this  case,  Theorem  6  states: 

N—l  N—l 

(fo.GSlVpOln]}  =  V  l(t>0,Jf[n']>  =  (v„,  £  * [»']}, 

n'= 0  n'= 0 

meaning  the  DC  component  of  every  entry  of  the  GST  is  the  DC  component  of  the 
sum  of  all  the  video  frames.  Such  lossy  behavior  is  accentuated  in  the  extreme  case 
where  U  =  I .  Here, 

iV-l 

GST jX[n]  =  x[n'}, 

n'= 0 

meaning  every  entry  of  the  GST  is  a  sum  of  all  the  entries  of  X. 

As  these  examples  illustrate,  the  idea  of  computing  a  GST  as  a  system  of  M 
DFTs  made  explicit  in  Theorem  6  is  a  useful  theoretical  concept.  However,  we  must 
point  out  that  it  offers  no  real  advantages  from  a  computational  point  of  view.  This 
method  of  computing  the  GST  requires  us  to  first  know  the  eigenvalues  and  eigen¬ 
vectors  of  U .  When  U  is  a  rotation  operator,  we  do  not  have  these  eigenvalues  and 
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eigenvectors  explicitly.  Moreover,  even  if  these  values  were  known,  such  as  in  the 
case  where  U  is  a  translation  operator,  the  lengthy  computation  time  of  this  method 
makes  it  undesirable. 

Indeed,  the  Theorem  6  method  for  finding  the  GST  requires  us  to  first  compute 
MN  inner  products  of  the  form  (vm,X\n'\)  in  HI,  where  HI  is  M-dimensional.  Thus, 
in  general,  when  the  v'ms  are  not  sparse,  we  expect  to  need  0(M2N )  operations  just 
to  compute  these  inner  products.  We  then  must  compute  M  FFTs  of  size  N  to 
find  the  values  (vm,  GSTj/(X)[n]),  requiring  0(MN  log  N)  operations  total.  We  then 
reconstruct  GST;yW[n]  according  to: 

M- 1 

GSTtW[n]  =  Y,  ( v m,  GST vX[n])vm,  (20) 

m= 0 

which  uses  an  additional  0(A42N)  operations,  for  a  total  complexity  of  0(M2N  + 
MN  log  N).  If  we  assume  U  has  sparse  eigenvectors,  this  value  decreases  to  0(MN  + 
MN  log  N). 

Meanwhile,  computing  the  GST  using  the  FGST  algorithm  of  Theorem  5  re¬ 
quires  0(N  log  N)  matrix- vector  products.  The  imrotate  command  in  Matlab  ap¬ 
plies  a  sparse  neighbor  interpolation  to  complete  each  matrix-vector  product,  each  re¬ 
quiring  O(M)  operations.  Hence,  the  total  computation  requirement  for  this  method 
is  O (MN  log  N).  Meanwhile,  computing  the  GST  directly  from  its  definition  requires 
0(N2)  matrix-vector  products,  and  is  therefore  the  slowest  method  overall.  Note  that 
although  the  methods  of  Theorems  5  and  6  both  require  O  (MN  log  N)  operations, 
the  FGST  is  less  complicated  and  does  not  require  the  eigenvalues  and  eigenvectors 
of  U  to  be  known. 

A  similar  comparison  can  be  made  when  U  is  a  circular  translation  operator, 
namely,  when  HI  =  £( Zjy),  U  =  T1,  and  M  =  N.  Here,  we  first  need  (vm,X\ n'])  for  all 
m,  n'  where  the  nm’s  are  the  Fourier  basis,  namely  the  eigenvectors  of  U  =  T.  For  any 
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fixed  n' ,  this  task  reduces  to  computing  the  DFT  of  X[n'].  Indeed,  since  the  number 
of  n! s  is  equal  to  N,  and  the  computational  cost  of  each  DFT  is  (9  (IV  log  IV). 

We  then  use  Theorem  6  to  compute  {vm,  GSTj/(AI)[n])  for  all  choices  of  m 
and  n,  incurring  an  additional  cost  of  N  DFTs  of  order  N.  Next,  we  reconstruct 
according  to  (20)  using  an  additional  N  DFTs  of  size  N.  Therefore,  in  the  case  of 
the  translation  operator,  the  total  number  of  operations  use  to  compute  the  GST 
according  to  Theorem  6  is  0(N2  log  IV). 

Meanwhile,  applying  the  FGST  algorithm  of  Theorem  5  to  the  circular  trans¬ 
lation  operator  requires  0(N  log  N)  matrix-vector  products,  each  of  which  can  be 
implemented  in  0(N )  operations  using  Matlab’s  circshift  command.  These  circ- 
shifts  only  involve  data  transfer  in  memory,  that  is,  no  multiplications.  Together, 
the  total  cost  of  computing  a  GST  using  Theorem  5  when  U  is  a  circular  translation 
operator  is  thus  0(N2  log IV). 

Although  the  methods  of  Theorems  5  and  6  are  comparable  in  terms  of  order 
of  operations,  the  FGST  is  less  complicated  and  involves  no  complex  arithmetic.  In 
conclusion,  we  see  that  in  every  case  in  which  we  wish  to  compute  the  GST  it  is 
better,  from  a  computational  point  of  view,  to  generalize  the  FFT  algorithm  in  the 
form  of  our  FGST  algorithm  (Theorem  5)  than  to  write  a  GST  as  a  system  of  FFTs 
in  the  spectral  domain  (Theorem  6).  In  fact,  in  the  applications  discussed  in  the 
next  chapter,  we  use  an  FGST  with  Matlab’s  imrotate  command  directly  without 
needing  to  have  the  eigenvalues  and  eigenvectors  of  our  rotation  operator. 
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VI.  Applications 

We  now  consider  several  examples  in  which  we  apply  the  FGST  algorithm  of  Theorem 
5  in  order  to  estimate  angular  frequency.  Consider  our  first  application  of  the  GST. 
a  pinwheel  rotating  at  a  constant  unknown  rate.  To  obtain  this  experimental  data, 
we  set  up  a  digital  camera  in  front  of  a  pinwheel  that  was  rotating  at  a  constant  rate 
with  air  flowing  from  a  nearby  fan. 

We  first  shot  a  series  of  still  frames  in  continuous  shooting  mode.  According 
to  Landau  in  [7],  every  finite  energy  signal  with  bandwidth  W  Hz  can  be  completely 
recovered  by  taking  samples  at  the  rate  of  2 W  per  second,  which  is  known  as  the 
Nyquist  rate.  However,  as  soon  as  the  GST  was  applied  to  this  image  sequence  it  was 
evident  by  every  GST  image  appearing  blurry  that  the  continuous  shooting  frame 
rate  did  not  meet  the  Nyquist  requirement. 

We  then  recorded  6  seconds  of  digital  video  and  converted  it  into  a  sequence  of 
192  image  frames  in  Matlab.  The  top  left  image  of  Figure  5  shows  an  example  of  one 
of  these  pinwheel  images.  After  computing  the  GST  of  this  sequence — computing  the 
GST  of  each  of  the  red,  green,  and  blue  channels  separately  an  then  combining  the 
result — the  frames  5-7  GST  images  and  frames  10-13  GST  images  are  examples  of 
blurry  GST  images.  Note  the  GST  images  from  frames  1-4  and  frames  13-192  are  not 
shown  since  they  look  similar  to  the  blurry  images  already  displayed.  These  blurry 
images  with  no  color  distinction  between  the  shades  indicate  the  rate  of  rotation  does 
not  match  that  particular  choice  of  n. 

In  the  frame  8  and  9  GST  images  however,  the  GST  is  are  clearly  sharper  and  has 
separated  colors.  It  makes  intuitive  sense  that  the  two  sharper  GST  images  correspond 
to  the  true  rate  of  rotation,  because  the  colors  have  been  added  coherently  on  top 
of  each  other.  During  this  video  sequence,  we  visually  estimated  that  the  pinwheel 
completed  8.5  revolutions  in  6  seconds.  Then,  the  GST  algorithm  corroborated  this 
as  the  total  number  of  revolutions  of  the  pinwheel,  as  seen  in  the  location  of  the  peak 
in  Figure  6. 
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Figure  5:  On  the  top  left,  the  pinwheel  image  shown  is  one  of  the  192-frame  pinwheel 
video  sequence.  The  frame  8  and  9  GST  images  correspond  to  the  true  rate 
of  rotation  of  the  pinwheel,  while  the  other  blurry  GST  images  do  not. 

For  our  second  GST  example,  we  obtained  a  100-frame  video  sequence  of  real 
data  from  a  200- Watt  thruster,  as  discussed  in  the  introduction.  This  thruster  was 
operating  at  80  Volts,  and  Liu  in  [8]  used  a  camera  to  capture  these  images  at  one 
million  frames  per  second  with  an  exposure  time  of  750  nanoseconds.  The  top  left 
image  in  Figure  7  is  an  example  of  plasma  being  ejected  from  the  thruster.  Each 
image  from  the  thruster  sequence  is  312  x  260  pixels,  and  in  light  of  Figure  4,  we  did 
not  attempt  to  compute  the  GST  directly.  On  a  standard  laptop  with  a  2.26  GHz 
Intel  Core  2  Duo  processor,  it  took  36  seconds  to  compute  the  GST  using  the  FGST 
algorithm. 

After  computing  the  GST  of  this  sequence  of  thruster  images,  the  frame  1-7 
GST  images  and  frame  9-14  GST  images  are  examples  of  blurry  GST  images.  Like 
the  pinwheel  example,  these  blurry  images  indicate  the  rate  of  rotation  does  not 
match  that  particular  choice  of  n.  While  challenging  to  verify  visually,  the  frame  8 
GST  image  corresponds  to  our  true  rate  of  plasma  rotation,  and  is  confirmed  by  the 
peak  of  the  norms  shown  in  Figure  8.  Interpreting  this  result  implies  the  plasma 
was  rotating  at  approximately  80  Kilohertz,  which  is  reasonable  given  the  parameters 
specified  in  [8]. 
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Figure  6:  The  P2  norm  and  the  £°°  norm  of  the  GST  of  the  pinwheel  video  sequence.  In 
both  norm  plots  there  is  a  peak  at  approximately  8.5,  corresponding  to  the 
number  of  revolutions  the  pinwheel  completed  in  6  seconds.  Note,  we  do  not 
show  the  graph  of  the  i1  norm  because  it  does  not  provide  any  information 
about  the  rate  of  rotation.  Indeed,  since  the  GST  is  computed  by  summing 
up  rotations  of  the  pinwheel  frames,  the  total  sums  as  calculated  by  the  ll 
norm  are  equal  in  all  GSTs  of  the  video  sequence. 


As  these  examples  have  shown,  the  GST  can  be  applied  to  real  data  and  provide 
an  accurate  estimate  of  the  rate  of  rotation.  While  not  shown  in  this  chapter,  the 
GST  can  also  be  applied  to  translation  problems  to  provide  an  accurate  estimate  of 
the  lateral  velocity.  In  fact,  the  appendix  includes  MATLAB  code  already  prepared 
with  circshift  commands,  so  that  is  can  be  applied  to  translation  problems.  Further 
guidance  is  also  included  in  the  appendix  about  how  to  change  the  code  appropriately 
for  rotation  problems. 
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Thruster 


Frame  1  GST  Frame  2  GST  Frame  3  GST  Frame  4  GST 
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**  O  O  O  C 


Frame  10  GST  Frame  11  GST  Frame  12  GST  Frame  13  GST  Frame  14  GST 

Figure  7:  On  the  top  left,  the  thruster  image  shown  is  one  of  the  100-frame  thruster 
sequence  collected  by  Liu  in  [8].  The  frame  8  GST  corresponds  to  the  true 
rate  of  rotation  of  the  plasma,  while  the  other  blurry  GST  images  do  not. 
The  frame  15-192  GSTs  are  not  shown  since  they  look  blurry  just  like  the 
ones  already  displayed. 


£2  Norm  £°°  Norm 

Figure  8:  The  £2  norm  and  the  £°°  norm  of  the  GST  of  the  thruster  video  sequence 
are  shown  above.  In  both  norm  graphs  there  is  a  peak  at  approximately  8, 
indicating  the  true  rate  of  rotation  of  the  plasma  is  80  Kilohertz. 
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Algorithm  for  Matlab  Code 

The  following  Matlab  code  shows  algorithms  for  the  direct  GST  and  the  FGST. 
The  direct  GST  should  not  be  implemented  since  it  is  not  as  computational  efficient 
as  the  FGST,  but  is  presented  below  for  comparison.  In  both  versions  of  the  code, 
X  is  a  video  sequence  (3D  data  cube),  and  must  be  defined  by  the  user  before  the 
function  can  be  called.  Both  code  versions  can  also  be  easily  generalized  to  the  case  of 
cyclic  translations  provided  the  imrotate  commands  are  replaced  by  uses  of  circshift. 

The  following  code  computes  the  GST  directly: 

N  =  size(X,3) ; 

Y  =  zeros(size(X)) ; 
for  n  =  0:N-1 
for  m  =  0:N-1 

Y(:,:,n+1)  =  Y( : , : ,n+l)+imrotate(X( : , : ,m+l) ,-l*m*n) ; 

end 

end 


This  second  set  of  code  computes  the  FGST  according  to  the  algorithm  of  Theorem  5 
and  is  implemented  by  only  using  a  small  amount  of  memory: 

°/„  Computing  the  prime  factorization  of  N  and  resulting  partial  products 
7.  according  to  Definition  4. 

CalP  =  factor(size(X,3)) ’ ; 

J  =  length (CalP) ; 

CalN  =  [  f lipud(cumprod(f lipud  (CalP)))  ;  1]; 

7„  Computing  the  FGST  in  bit  reversal  order  according  to  Theorem  5. 
CumProd  =  [1;  cumprod(CalP)] ; 

B  =  0; 

for  j  =  0 : J— 1 ; 
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%  Computing  the  bit  reversal  permutation  (8)  according  to  (9). 

B  =  kron(B,  ones(CalP(j+l) , 1) ) 

+CumProd( j+1) *kron(ones (CumProd( j+1) , 1) , (0 : CalP( j+l)-l) 5 ) ; 

for  k  =  0:CalN(j+2)-l; 

for  1  =  0 : CalN(l) /CalN (j+1)  -1; 

Indices  =  CalN(j+2)*(CalP(j+l)*l+(0:CalP(j+l)  -l))+k; 

XI  =  X ( : , : , Indices+1) ; 
for  p  =  0 : CalP (j+1)  -1; 

X( : , : , Indices (p+l)+l)=zeros (size (X ( : , : , Indices (p+l)+l))) ; 
for  q  =  0 : CalP( j+l)-l ; 

X( : , : ,Indices(p+l)+l)  =  X( Indices (p+l)+l) 
+imrotate(Xl ( : , : ,q+l) ,-CalN(j+2)*q*B(CalP(j+l)*l+p+l)) 
end 

end 

end 

end 

end 

7.  Permuting  the  bit-reversed  FGST  into  standard  order, 
for  n=0:CalN(l)-l; 
index=  find(B==n); 

XI  =  X( : , : , index) ; 

X(:,:, index)  =  X(:,:,n+1); 

X(: , : ,n+l)  =  XI; 

B(index)  =  B(n+1); 

B(n+1)  =  n; 

end 
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